THE HEMODYNAMIC RESPONSE FUNCTION VARIES ACROSS ANATOMICAL LOCATION AND PATHOLOGY IN THE EPILEPTIC BRAIN (SUPPLEMENTARY: CODE)

Zhengchen Cai¹, Nicolás Von Ellenrieder¹, Thaera Arafat¹, Hui Ming Kho², Gang Chen³, Andreas Koupparis⁴, Chifaou Abdallah¹, Roy Dudley⁵, Dang Khoa Nguyen⁶, Jeff Hall¹, Francois Dubeau¹, Jean Gotman¹, Boris Bernhardt¹

1. The Neuro (Montreal Neurological Institute-Hospital), McGill University, Montreal, Canada.

2. Department of Neurosurgery, Osaka University Graduate School of Medicine, Suita, Japan.

3. Scientific and Statistical Computing Core, National Institute of Mental Health, Bethesda, USA.

4. The Cyprus Institute of Neurology and Genetics, Nicosia, Cyprus.

5. Montreal Children’s Hospital, McGill University, Montreal, Canada.

6. Centre Hospitalier de l’Université de Montréal, Montréal, Canada.

Preface

This document contains the code used to produce the results in the paper. Each results section contains the corresponding computation code chunk. All code chunks and printed messages are folded by default. They can be unfolded to see the details. This provides a more coherent way of presenting results and code together.

Library

Library used in the analysis are listed below.

library(tidyverse)
library(ggplot2)
library(stringr)
library(ggpubr)
library(scales)
library(ragg)
library(readxl)
library(gt)
library(ggstatsplot)
library(brms)
library(rstan)
library(tidybayes)
library(webr)
library(DT)
library(dagitty)
library(ggdag)
library(kableExtra)
library(moonBook)
library(ggforce)
library(ggsci)
library(marginaleffects)
library(ggrepel)
library(brainspriteR)
library(readr)
library(nnet)
library(neuRosim)
library(lme4)
library(mclogit)
library(ggseg)
library(ggsegDesterieux)
library(ggsegSchaefer)
library(ggridges)
library(tidymodels)  
library(themis)      
library(xgboost)     
library(ranger)      
library(tidyr)
library(vip)

Data wrangling

Load data and prepare variable values.

derivatives_dir <- "/Volumes/django/EP_EEGfMRI_deconvolution/derivatives"

output_dir <- "/Volumes/django/EP_EEGfMRI_deconvolution/derivatives/results"

# The neuRosim::canonicalHRF() implements the standard double‑gamma haemodynamic response function as originally characterized in event‑related fMRI by Friston et al. (1998) and refined by Glover (1999)

# Parameter Description Default
# a1    Delay of the main (positive) response   6s
# a2    Delay of the undershoot 12s
# b1    Dispersion (width) of the main response 0.9
# b2    Dispersion of the undershoot    0.9
# c Scale (height) of the undershoot relative to the main response  0.35

hrfs <- read_csv('./hrflib18.csv', show_col_types = FALSE) 

hrf_pcs <- read_csv('./HRFbasis_500Hz.csv', show_col_types = FALSE) %>% 
  rename(PC1 = pc1, PC2 = pc2, PC3 = pc3) %>% 
  arrange(time) %>% 
  pivot_longer(
    cols      = -time,        
    names_to  = "pc",
    values_to = "value"
  )

# default 
hrf_canonical <- canonicalHRF(hrfs$hrf_time)    

long_hrf <- hrfs %>% 
  arrange(hrf_time) %>% 
  pivot_longer(
    cols      = -hrf_time,        
    names_to  = "hrf",
    values_to = "value"
  ) %>% 
  group_by(hrf) %>% 
  mutate(
    corr = cor(value, hrf_canonical),
    rmse = sqrt(mean((value - hrf_canonical)^2))
  ) %>% 
  ungroup() %>% 
  mutate(
    hrf_group = case_when(
      hrf == "hrf1" ~ "group2",
      hrf == "hrf2" ~ "group2",
      hrf == "hrf3" ~ "group1",
      hrf == "hrf4" ~ "group2",
      hrf == "hrf5" ~ "group1",
      hrf == "hrf6" ~ "group2",
      hrf == "hrf7" ~ "group1",
      hrf == "hrf8" ~ "group1",
      hrf == "hrf9" ~ "group1",
      hrf == "hrf10" ~ "group3",
      hrf == "hrf11" ~ "group3",
      hrf == "hrf12" ~ "group3",
      hrf == "hrf13" ~ "group4",
      hrf == "hrf14" ~ "group3",
      hrf == "hrf15" ~ "group4",
      hrf == "hrf16" ~ "group4",
      hrf == "hrf17" ~ "group4",
      hrf == "hrf18" ~ "group4",
    ) ) %>% 
  mutate(panel = hrf_group) %>% 
  group_by(hrf_group) %>% 
  mutate(corr_group = mean(corr),
         rmse_group = mean(rmse)) %>% 
  ungroup()

cor_canonical <- long_hrf %>% 
 distinct(hrf_group, corr_group, rmse_group) %>% 
 arrange(hrf_group)


ylim_hrfs <- range(long_hrf$value) * 1.05

panel_labels <- c(
  group1 = "group1: no undershoot",
  group2 = "group2: low undershoot",
  group3 = "group3: moderate and early undershoot",
  group4 = "group4: large and later undershoot"
)

long_hrf %>% 
  ggplot(aes(x = hrf_time, y = value, group = hrf, color = rmse)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(color = "white", size = 1.6) +
  geom_line(size = 0.8, alpha = 1) +
    scale_color_gradientn(
    name   = "RMSE to the canonical HRF",
    colours = c(
      "#fee0d2",  
      "#fcbba1",
      "#fc9272",
      "#ef3b2c",
      "#cb181d",
      "#a50f15",
      "#67000d"   
    ),
    guide   = guide_colorbar(
      title.position = "top",
      title.hjust    = 0.5,
      barwidth       = 15,
      barheight      = 1
    )
  ) +
  scale_y_continuous(
    breaks = 0,
    limits = ylim_hrfs
  ) +
  facet_wrap(
    ~ panel,
    ncol     = 2,
    labeller = labeller(panel = panel_labels)
  ) +
  labs(
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    strip.text.x        = element_text(hjust = 0, size = 18),
    axis.title.x        = element_text(size = 18),
    axis.text.x         = element_text(size = 18),  
    axis.text.y         = element_text(size = 14),
    panel.grid.major.x  = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y  = element_blank(),
    panel.grid.minor    = element_blank(),
    panel.spacing       = unit(1, "lines"),
    legend.position     = "bottom",
    legend.direction    = "horizontal",
    legend.title        = element_text(size = 16),
    legend.text         = element_text(size = 16)    
  ) -> p 

p_file1 <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_file1, width = 10, height = 16/2, units = "in", res = 300)
suppressWarnings(print(p))
invisible(dev.off())


anno_df <- hrf_pcs %>%
  distinct(pc) %>%
  mutate(x = 10, y = 0)

hrf_pcs %>% 
  ggplot(aes(x = time, y = value, group = pc)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(size = 0.8, alpha = 1) + 
  scale_y_continuous(breaks = 0) +
  facet_wrap(~ pc, ncol = 1) +
  geom_text(
    data = anno_df,
    aes(x = x, y = y, label = pc),
    size = 6,
    hjust = 0
  ) +
  labs(
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    strip.text           = element_blank(),           
    strip.background     = element_blank(),
    axis.title.x         = element_text(size = 20),
    axis.text.x          = element_text(size = 20),
    axis.text.y          = element_text(size = 20),
    panel.grid.major.x   = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y   = element_blank(),
    panel.grid.minor     = element_blank(),
    panel.spacing        = unit(0.2, "lines"),         
    legend.position      = "none"
  ) -> p1
  
long_hrf %>% 
  ggplot(aes(x = hrf_time, y = value, group = hrf, color = hrf_group)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(color = "white", size = 1.6) +
  geom_line(size = 0.8, alpha = 0.6) +
  scale_color_manual(values = c("#00468BFF", "#00468BFF", "#00468BFF","#00468BFF")) + 
  scale_y_continuous(
    breaks = 0,
    limits = ylim_hrfs
  ) +
  labs(
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    strip.text.x        = element_text(hjust = 0, size = 20),
    axis.title.x        = element_text(size = 20),
    axis.text.x         = element_text(size = 20),   
    axis.text.y         = element_text(size = 20),
    panel.grid.major.x  = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y  = element_blank(),
    panel.grid.minor    = element_blank(),
    panel.spacing       = unit(1, "lines"),
    legend.position     = "none",
    legend.direction    = "horizontal",
    legend.title        = element_blank(),
    legend.text         = element_text(size = 16)    
  ) -> p2

p <- ggarrange(p1, p2, nrow = 1, ncol = 2, widths = c(1, 1))

p_file2 <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_file2, width = 10, height = 4.2, units = "in", res = 300)
suppressWarnings(print(p))
invisible(dev.off())

# cavity_tibble contains all voxel xyz inside the cavity per subj per surgery 

cavity_tibble <- read_csv('./cavity_tibble.csv', show_col_types = FALSE) %>% 
    rename(postop_date = surgery) %>% 
    mutate(postop_date = as.Date(postop_date, format = "%d_%m_%Y")) %>%
    group_by(subj, voxel) %>%
    slice_min(postop_date, n = 1, with_ties = FALSE) %>%
    ungroup() %>% 
    mutate(postop_date = format(postop_date, "%d_%m_%Y"))

# hrflib contains all HRF information per voxel

hrflib <- read_csv("./table_hrflib18_withsign.csv", show_col_types = FALSE) %>% 
    mutate(type_dur = ifelse(str_detect(type_name, "-dur"), "long-lasting", "spike")) %>% 
    mutate(subj = str_replace(subj, "sub_", "sub-")) 

fs_label <- read_csv("./afni_fs_aparc+aseg_2009_goodnames.csv", show_col_types = FALSE) %>% 
    rename(fs_label = AFNI_label) %>% 
    mutate(lobe = str_c(hemisphere, " ", lobe))

# ep_info contains clincal information with surgery date

ep_info <- read_csv("./ep_info.csv", show_col_types = FALSE) %>% 
  mutate(age_scan = round(age - as.numeric(difftime(surgery_date, scan_date, units = "days")) / 365),
         duration_of_epilepsy_scan = round(duration_of_epilepsy - as.numeric(difftime(surgery_date, scan_date, units = "days")) / 365),
         surgery_date = format(surgery_date, "%d_%m_%Y"),
         postop_date = format(postop_date, "%d_%m_%Y")) %>% 
  rename(subj = participant_id,
         surgery = surgery_date) %>% 
  group_by(subj) %>% 
  mutate(age_scan = round(mean(age_scan)), 
         duration_of_epilepsy_scan = round(mean(duration_of_epilepsy_scan))) %>% 
  ungroup() %>% 
  select(subj, gender, age_scan, duration_of_epilepsy_scan, pathology, surgery, postop_date, MRIlesion)

# this step match the clinical info particular to patient and surgery (by date)
cavity_tibble <- cavity_tibble %>% 
  left_join(ep_info, by = c("subj", "postop_date")) %>% 
  select(-c("postop_date"))


# add voxel label 
hrflib <- hrflib %>% 
    mutate(
    hrf_group = case_when(
      hrf_id == "1" ~ "2",
      hrf_id == "2" ~ "2",
      hrf_id == "3" ~ "1",
      hrf_id == "4" ~ "2",
      hrf_id == "5" ~ "1",
      hrf_id == "6" ~ "2",
      hrf_id == "7" ~ "1",
      hrf_id == "8" ~ "1",
      hrf_id == "9" ~ "1",
      hrf_id == "10" ~ "3",
      hrf_id == "11" ~ "3",
      hrf_id == "12" ~ "3",
      hrf_id == "13" ~ "4",
      hrf_id == "14" ~ "3",
      hrf_id == "15" ~ "4",
      hrf_id == "16" ~ "4",
      hrf_id == "17" ~ "4",
      hrf_id == "18" ~ "4",
    ) ) %>% 
   left_join(fs_label %>% select(fs_label, labelname, lobe, tissue), by = "fs_label") %>% 
   left_join(cavity_tibble %>% distinct(subj, gender, age_scan), by = "subj")

# add pathology
hrflib_clean <- hrflib %>% 
   left_join(cavity_tibble %>% distinct(subj, voxel, engel, pathology, surgery, MRIlesion), by = c("subj", "voxel")) %>% 
   mutate(pathology_complex = pathology) %>% 
   mutate(pathology = case_when(
      pathology_complex == "Cortical malformation" ~ "cortical malformation",
      pathology_complex == "non lesional" ~ "normal",
      pathology_complex == "HS" ~ "HS|dual",
      pathology_complex == "dual pathology" ~ "HS|dual",
      pathology_complex == "Gliosis" ~ "gliosis",
      pathology_complex == "encephalitis" ~ "enceph|tumor|vascular",
      pathology_complex == "Vascular malformation" ~ "enceph|tumor|vascular",
      pathology_complex == "complex pathology" ~ "complex",
      pathology_complex == "Tumor" ~ "enceph|tumor|vascular",
      TRUE ~ NA)) %>% 
   rename(age = age_scan)  %>% 
   filter(abs(r) >= 0.7,
         !is.na(lobe), 
         tissue == "tiss__gm") %>% 
   filter(!fs_label %in% c("6", "13", "26")) %>% 
   filter(
    abs(posPeakAmp) < 2.5,
    posPeakTime > 3 & posPeakTime < 13,
    abs(negPeakAmp) < 2.5,
  ) 


df_all <- hrflib_clean %>%
  group_by(subj, type, fs_label) %>%
  mutate(
    pathology = {
      non_na_path <- pathology[!is.na(pathology)]
      if (length(non_na_path) == 0) {
        NA
      } else {
        most_common <- names(sort(table(non_na_path), decreasing = TRUE))[1]
        most_common
      }
    }
  ) %>%
  ungroup() %>% 
  filter(is.na(pathology) | pathology == "normal") %>%  
  group_by(fs_label) %>%
  mutate(voxelinparcel = n()) %>%
  ungroup() %>% 
  filter(voxelinparcel >= 5) %>% 
  mutate(sign = ifelse(sign == 1, "positive", "negative")) %>% 
  mutate(gender = factor(gender),
         fs_label = factor(fs_label),
         labelname = factor(labelname),
         subj = factor(subj),
         lobe = factor(lobe),
         sign = factor(sign)) %>% 
  mutate(hrf_id = factor(hrf_id),
         hrf_group = factor(hrf_group, levels = c("1", "2", "3", "4"))) %>%
  select(subj, hrf_id, hrf_group, age, gender, type_dur, lobe, fs_label, labelname,
         posPeakAmp, posPeakTime, posPeakFWHM, negPeakAmp, negPeakTime, negPeakFWHM, sign)


df <- hrflib_clean %>%
  filter(!is.na(pathology)) %>% 
  group_by(fs_label) %>%
  mutate(voxelinparcel = n()) %>%
  ungroup() %>% 
  filter(voxelinparcel >= 5) %>% 
  group_by(pathology_complex) %>%
  mutate(voxelinparcel = n()) %>%
  ungroup() %>% 
  filter(voxelinparcel >= 5) %>% 
  mutate(sign = ifelse(sign == 1, "positive", "negative")) %>% 
  mutate(pathology = ifelse(pathology == "enceph|tumor|vascular", "encephalitis", as.character(pathology))) %>% 
  mutate(gender = factor(gender),
         fs_label = factor(fs_label),
         labelname = factor(labelname),
         subj = factor(subj),
         lobe = factor(lobe),
         pathology = factor(pathology),
         sign = factor(sign)) %>% 
  mutate(hrf_id = factor(hrf_id),
         hrf_group = factor(hrf_group, levels = c("1", "2", "3", "4"))) %>% 
  select(subj, hrf_id, hrf_group, age, gender, type_dur, lobe, fs_label, pathology, labelname,
         posPeakAmp, posPeakTime, posPeakFWHM, negPeakAmp, negPeakTime, negPeakFWHM, sign)


p_file_posPeakAmp <-
  plot_histrogram(df_all, "posPeakAmp", "BOLD % change", title = bquote("distribution of HRF peak amplitude"), NULL)

p_file_posPeakTime <-
  plot_histrogram(df_all, "posPeakTime", "Time(s)", title = bquote("distribution of HRF time to peak"), NULL)

p_file_posPeakFWHM <-
  plot_histrogram(df_all, "posPeakFWHM", "FWHM(s)", title = bquote("distribution of HRF peak FWHM"), NULL)

p_file_negPeakAmp <-
  plot_histrogram(df_all, "negPeakAmp", "BOLD % change", title = bquote("distribution of HRF undershoot amplitude"), NULL)

p_file_negPeakTime <-
  plot_histrogram(df_all %>% filter(negPeakTime>7), "negPeakTime", "Time(s)", title = bquote("distribution of HRF time to undershoot"), NULL)

p_file_negPeakFWHM <-
  plot_histrogram(df_all, "negPeakFWHM", "FWHM(s)", title = bquote("distribution of HRF undershoot FWHM"), NULL)

HRF lib

hrf_group	RMSE to the canonical HRF
group1	0.55
group2	0.37
group3	0.23
group4	0.56

Summary of the variables

Bayesian analysis

Within less pathological regions

HRF features

Peak amplitude

mean_posPeakAmp = mean(df_all$posPeakAmp)
sd_posPeakAmp = sd(df_all$posPeakAmp)

fit2 <- brm(
  data = df_all %>% mutate(posPeakAmp = (posPeakAmp - mean_posPeakAmp)/sd_posPeakAmp),
  family = gaussian(),
  formula = posPeakAmp ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakAmp_nopathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakAmp_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_posPeakAmp) + mean_posPeakAmp)

roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

colorbar_lim <- c(0.04, 0.56)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50",
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "a | HRF peak amplitude") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "BOLD % change",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50",
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1


proisc1a <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1a))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Time to peak

mean_posPeakTime = mean(df_all$posPeakTime)
sd_posPeakTime = sd(df_all$posPeakTime)

fit2 <- brm(
  data = df_all %>% mutate(posPeakTime = (posPeakTime - mean_posPeakTime)/sd_posPeakTime),
  family = gaussian(),
  formula = posPeakTime ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakTime_nopathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakTime_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_posPeakTime) + mean_posPeakTime)


roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

colorbar_lim <- c(floor(min(roi_temp$estimate)*10)/10, ceiling(max(roi_temp$estimate)*10)/10)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "b | HRF time to peak") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "Time(s)",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1


proisc1b <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1b))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Peak width

mean_posPeakFWHM = mean(df_all$posPeakFWHM)
sd_posPeakFWHM = sd(df_all$posPeakFWHM)

fit2 <- brm(
  data = df_all %>% mutate(posPeakFWHM = (posPeakFWHM - mean_posPeakFWHM)/sd_posPeakFWHM),
  family = gaussian(),
  formula = posPeakFWHM ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakFWHM_nopathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakFWHM_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_posPeakFWHM) + mean_posPeakFWHM)


roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()

colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

colorbar_lim <- c(floor(min(roi_temp$estimate)*10)/10, ceiling(max(roi_temp$estimate)*10)/10)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "c | HRF peak FWHM") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "FWHM(s)",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1

subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    limits = colorbar_lim,
    na.value = "grey50", 
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1


proisc1c <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1c))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Undershoot amplitude

mean_negPeakAmp = mean(df_all$negPeakAmp)
sd_negPeakAmp = sd(df_all$negPeakAmp)

fit2 <- brm(
  data = df_all %>% mutate(negPeakAmp = (negPeakAmp - mean_negPeakAmp)/sd_negPeakAmp),
  family = gaussian(),
  formula = negPeakAmp ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakAmp_nopathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakAmp_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_negPeakAmp) + mean_negPeakAmp)


roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()


colorbar_col <- adjustcolor(
  rev(c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")),
  alpha.f = 1)

colorbar_lim <- c(-0.56, -0.04)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "d | HRF undershoot amplitude") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "BOLD % change",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1


proisc1d <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1d))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Time to undershoot

mean_negPeakTime = mean(df_all$negPeakTime)
sd_negPeakTime = sd(df_all$negPeakTime)

fit2 <- brm(
  data = df_all %>% mutate(negPeakTime = (negPeakTime - mean_negPeakTime)/sd_negPeakTime),
  family = gaussian(),
  formula = negPeakTime ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakTime_nopathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakTime_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_negPeakTime) + mean_negPeakTime)


roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

colorbar_lim <- range(roi_temp$estimate)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "e | HRF time to undershoot") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "Time(s)",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1

proisc1e <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1e))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Undershoot width

mean_negPeakFWHM = mean(df_all$negPeakFWHM)
sd_negPeakFWHM = sd(df_all$negPeakFWHM)

fit2 <- brm(
  data = df_all %>% mutate(negPeakFWHM = (negPeakFWHM - mean_negPeakFWHM)/sd_negPeakFWHM),
  family = gaussian(),
  formula = negPeakFWHM ~ age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakFWHM_nopathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakFWHM_nopathology.rds")

mfx <- avg_predictions(fit2, by = "fs_label") %>% 
  mutate(estimate = (estimate * sd_negPeakFWHM) + mean_negPeakFWHM)


roi_temp <- mfx %>% group_by(fs_label) %>%
  select(fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname) %>% 
  ungroup()


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

colorbar_lim <- c(floor(min(roi_temp$estimate)*10)/10, ceiling(max(roi_temp$estimate)*10)/10)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    na.value = "grey50", 
    limits = colorbar_lim,
    breaks = pretty_breaks(5)) +
  labs(subtitle = "f | HRF undershoot FWHM") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "FWHM(s)",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> proi1

subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn( 
    colours = colorbar_col, 
    limits = colorbar_lim,
    na.value = "grey50", 
    breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),          
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 30, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "top", 
                                title.hjust = 0.5,
                                label.position = "bottom")) -> psc1

proisc1f <- ggarrange(proi1, psc1, nrow = 1, ncol = 2, common.legend = T, 
                      legend="bottom", widths = c(6, 1))


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 4.2, units = "in", res = figdpi)
suppressWarnings(print(proisc1f))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

HRF shape

fit1 <- brm(
  data = df_all, family = categorical(link = "logit", refcat = "1"),
  formula = hrf_group ~ 0 + age + gender + type_dur + sign + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 3), class = b, dpar = mu2),
    prior(normal(0, 3), class = b, dpar = mu3),
    prior(normal(0, 3), class = b, dpar = mu4),
    prior(exponential(1), class = sd, dpar = mu2),
    prior(exponential(1), class = sd, dpar = mu3),
    prior(exponential(1), class = sd, dpar = mu4)
  ),
  control = list(adapt_delta = 0.9),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221)

saveRDS(fit1, "fit_brms_model_4groups_nopathology.rds")

fit1 <- readRDS("fit_brms_model_4groups_nopathology.rds")

mfx <- avg_predictions(fit1, by = "fs_label")


roi_temp <- mfx %>% group_by(fs_label) %>%
  slice_max(order_by = estimate, n = 1, with_ties = FALSE) %>%
  ungroup() %>% 
  select(group, fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname)


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  labs(subtitle = "b | posterior probability of selected HRF per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 24),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> proi1

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  labs(subtitle = "a | most probable HRF group per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),  
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 26, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF groups",
    title.position = "top",
    title.hjust = 0.5,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> proi2


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> psc1


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF group",
    title.position = "left",
    title.vjust = 1,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> psc2


proisc2a <- ggarrange(proi2, psc2, nrow = 1, ncol = 2, common.legend = T, 
                     legend="bottom", widths = c(6, 1))


panel_labels_text <- cor_canonical %>% 
  mutate(rmse = round(rmse_group, 3)) %>% 
  pull(rmse)
  

panel_labels <- c(
  group1 = paste0("RMSE to canonical HRF: ", panel_labels_text[1]),
  group2 = paste0("RMSE to canonical HRF: ", panel_labels_text[2]),
  group3 = paste0("RMSE to canonical HRF: ", panel_labels_text[3]),
  group4 = paste0("RMSE to canonical HRF: ", panel_labels_text[4])
)

long_hrf %>% 
  ggplot(aes(x = hrf_time, y = value, group = hrf, color = hrf_group)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(color = "white", size = 1.6) +
  geom_line(size = 0.8, alpha = 1) +
  scale_color_manual(values = c("#00468BFF", "#ED0000FF", "#42B540FF",  "#0099B4FF")) + 
  scale_y_continuous(
    breaks = 0,
    limits = ylim_hrfs
  ) +
  facet_wrap(
    ~ panel,
    ncol     = 4,
    labeller = labeller(panel = panel_labels)
  ) +
  labs(
    subtitle = "b | HRF groups",
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    strip.text = ggtext::element_markdown(),
    strip.text.x        = element_text(hjust = 0, size = 17),
    axis.title.x        = element_text(size = 24),
    axis.text.x         = element_text(size = 20),   
    axis.text.y         = element_text(size = 20),
    panel.grid.major.x  = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y  = element_blank(),
    panel.grid.minor    = element_blank(),
    panel.spacing       = unit(1, "lines"),
    legend.position     = "none",
    legend.direction    = "horizontal",
    legend.title        = element_text(size = 20),
    legend.text         = element_text(size = 20)    
  ) -> proisc2b 

proisc2 <- ggarrange(proisc2a, proisc2b, nrow = 2, ncol = 1, heights = c(1, 1), align = "hv")

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 2*4.5, units = "in", res = figdpi)
suppressWarnings(print(proisc2))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

mfx <- avg_predictions(fit1, by = c("fs_label", "sign"))


roi_temp <- mfx %>% filter(sign == "positive") %>% group_by(fs_label) %>%
  slice_max(order_by = estimate, n = 1, with_ties = FALSE) %>%
  ungroup() %>% 
  select(group, fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname)


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  labs(subtitle = "b | posterior probability of selected HRF per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 24),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> proi1

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  labs(subtitle = "a | most probable HRF group per ROI (activation)") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),  
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 26, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF groups",
    title.position = "top",
    title.hjust = 0.5,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> proi2


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> psc1


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF group",
    title.position = "left",
    title.vjust = 1,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> psc2


proisc2a <- ggarrange(proi2, psc2, nrow = 1, ncol = 2, common.legend = T, 
                      legend="bottom", widths = c(6, 1))

panel_labels_text <- cor_canonical %>% 
  mutate(rmse = round(rmse_group, 3)) %>% 
  pull(rmse)


panel_labels <- c(
  group1 = paste0("RMSE to canonical HRF: ", panel_labels_text[1]),
  group2 = paste0("RMSE to canonical HRF: ", panel_labels_text[2]),
  group3 = paste0("RMSE to canonical HRF: ", panel_labels_text[3]),
  group4 = paste0("RMSE to canonical HRF: ", panel_labels_text[4])
)

long_hrf %>% 
  ggplot(aes(x = hrf_time, y = value, group = hrf, color = hrf_group)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(color = "white", size = 1.6) +
  geom_line(size = 0.8, alpha = 1) +
  scale_color_manual(values = c("#00468BFF", "#ED0000FF", "#42B540FF",  "#0099B4FF")) + 
  scale_y_continuous(
    breaks = 0,
    limits = ylim_hrfs
  ) +
  facet_wrap(
    ~ panel,
    ncol     = 4,
    labeller = labeller(panel = panel_labels)
  ) +
  labs(
    subtitle = "b | HRF groups",
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    strip.text = ggtext::element_markdown(),
    strip.text.x        = element_text(hjust = 0, size = 17),
    axis.title.x        = element_text(size = 24),
    axis.text.x         = element_text(size = 20),   
    axis.text.y         = element_text(size = 20),
    panel.grid.major.x  = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y  = element_blank(),
    panel.grid.minor    = element_blank(),
    panel.spacing       = unit(1, "lines"),
    legend.position     = "none",
    legend.direction    = "horizontal",
    legend.title        = element_text(size = 20),
    legend.text         = element_text(size = 20)    
  ) -> proisc2b 

proisc2 <- ggarrange(proisc2a, proisc2b, nrow = 2, ncol = 1, heights = c(1, 1), align = "hv")

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 2*4.5, units = "in", res = figdpi)
suppressWarnings(print(proisc2))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

roi_temp <- mfx %>% filter(sign == "negative") %>% group_by(fs_label) %>%
  slice_max(order_by = estimate, n = 1, with_ties = FALSE) %>%
  ungroup() %>% 
  select(group, fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname)


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  labs(subtitle = "b | posterior probability of selected HRF per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 24),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> proi1

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  labs(subtitle = "a | most probable HRF group per ROI (deactivation)") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),  
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 26, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF groups",
    title.position = "top",
    title.hjust = 0.5,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> proi2


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> psc1


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF group",
    title.position = "left",
    title.vjust = 1,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> psc2


proisc2a <- ggarrange(proi2, psc2, nrow = 1, ncol = 2, common.legend = T, 
                      legend="bottom", widths = c(6, 1))


panel_labels_text <- cor_canonical %>% 
  mutate(rmse = round(rmse_group, 3)) %>% 
  pull(rmse)


panel_labels <- c(
  group1 = paste0("RMSE to canonical HRF: ", panel_labels_text[1]),
  group2 = paste0("RMSE to canonical HRF: ", panel_labels_text[2]),
  group3 = paste0("RMSE to canonical HRF: ", panel_labels_text[3]),
  group4 = paste0("RMSE to canonical HRF: ", panel_labels_text[4])
)

long_hrf %>% 
  ggplot(aes(x = hrf_time, y = value, group = hrf, color = hrf_group)) +
  geom_hline(yintercept = 0, color = "grey", size = 0.2) +
  geom_line(color = "white", size = 1.6) +
  geom_line(size = 0.8, alpha = 1) +
  scale_color_manual(values = c("#00468BFF", "#ED0000FF", "#42B540FF",  "#0099B4FF")) + 
  scale_y_continuous(
    breaks = 0,
    limits = ylim_hrfs
  ) +
  facet_wrap(
    ~ panel,
    ncol     = 4,
    labeller = labeller(panel = panel_labels)
  ) +
  labs(
    subtitle = "b | HRF groups",
    x = "Time (s)",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    strip.text = ggtext::element_markdown(),
    strip.text.x        = element_text(hjust = 0, size = 17),
    axis.title.x        = element_text(size = 24),
    axis.text.x         = element_text(size = 20),   
    axis.text.y         = element_text(size = 20),
    panel.grid.major.x  = element_line(color = "grey80", size = 0.2),
    panel.grid.major.y  = element_blank(),
    panel.grid.minor    = element_blank(),
    panel.spacing       = unit(1, "lines"),
    legend.position     = "none",
    legend.direction    = "horizontal",
    legend.title        = element_text(size = 20),
    legend.text         = element_text(size = 20)   
  ) -> proisc2b 

proisc2 <- ggarrange(proisc2a, proisc2b, nrow = 2, ncol = 1, heights = c(1, 1), align = "hv")

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 2*4.5, units = "in", res = figdpi)
suppressWarnings(print(proisc2))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

mfx <- avg_predictions(fit1, by = "sign") %>% 
  posterior_draws() %>% 
  select(sign, draw, drawid, group) %>% 
  mutate(sign = str_replace_all(sign, "negative", "deactivation")) %>% 
  mutate(sign = str_replace_all(sign, "positive", "activation"))

group_colors <- c(
  "1" = "#00468BFF",
  "2" = "#ED0000FF",
  "3" = "#42B540FF",
  "4" = "#0099B4FF"
)

mfx <- mfx %>%
  mutate(group = factor(group),
         group_colored = paste0(
           "<span style='color:", group_colors[as.character(group)], "'>",
           as.character(group), "</span>"
         ))

mfx$group_colored <- factor(mfx$group_colored, levels = unique(mfx$group_colored))

mfx %>%
  ggplot(aes(x = draw, y = group_colored, fill = after_stat(x > 1/4))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("gray80", "skyblue")) +
  geom_vline(xintercept = 0.25, linetype = "dashed") +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3, 0.4)) +
  facet_wrap(~sign, nrow = 2, ncol = 3, scales = "free_x") +  
  labs(subtitle = "c | marginal effect of sign on HRF group",
       x = "Probability",
       y = "HRF group") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 24),
    axis.text.x = element_text(size = 20),   
    axis.text.y = ggtext::element_markdown(size = 20, margin = margin(r = 20))
  ) -> p

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 18/1.6, height = 8/1.6, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Summary of the distribution

# A tibble: 8 × 8
# Groups:   sign, group [8]
  sign         group Parameter Median    CI CI_low CI_high    pd
  <chr>        <fct> <chr>      <dbl> <dbl>  <dbl>   <dbl> <dbl>
1 activation   1     Posterior 0.140   0.95 0.135    0.146     1
2 activation   2     Posterior 0.417   0.95 0.409    0.426     1
3 activation   3     Posterior 0.144   0.95 0.138    0.150     1
4 activation   4     Posterior 0.298   0.95 0.291    0.306     1
5 deactivation 1     Posterior 0.285   0.95 0.275    0.296     1
6 deactivation 2     Posterior 0.306   0.95 0.296    0.318     1
7 deactivation 3     Posterior 0.312   0.95 0.302    0.321     1
8 deactivation 4     Posterior 0.0964  0.95 0.0893   0.104     1

mfx <- fit1 %>% 
  avg_slopes(variables = "age") %>% 
  posterior_draws() %>% 
  select(draw, drawid, group) 

mfx %>%
  mutate(draw = 500 * draw) %>% 
  ggplot(aes(x = draw, y = group, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), fatten_point = 1.8,
               point_interval = "median_hdi", normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("gray80", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(title = "marginal effect of age on HRF group",
       x = "probability in %",
       y = "HRF group",
       subtitle = "\u2206Pr(HRF group change) per 5y \u2206 in age") +
  theme(
    legend.position = "none",
    plot.margin = margin(5.5, 50, 5.5, 5.5, "pt"),
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(8, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = element_text(size = 32),
    text = element_text(size = 36, face = "plain", color = "black"),
    axis.text = element_text(color = "black", size = 36),
    axis.text.y = element_text(margin = margin(r = 20)),
    plot.title.position = "plot",
    plot.title = element_text(size = 46),
    plot.subtitle = element_text(size = 40),
    plot.caption = element_text(size = 40)) -> p

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 18, height = 12, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Summary of the distribution

# A tibble: 4 × 7
# Groups:   group [4]
  group Parameter Median    CI CI_low CI_high    pd
  <fct> <chr>      <dbl> <dbl>  <dbl>   <dbl> <dbl>
1 1     Posterior  -3.16  0.95  -3.44  -2.90      1
2 2     Posterior  -1.34  0.95  -1.71  -0.987     1
3 3     Posterior   2.91  0.95   2.61   3.20      1
4 4     Posterior   1.59  0.95   1.29   1.91      1

mfx <- avg_predictions(fit1, by = "gender") %>% 
  posterior_draws() %>% 
  select(gender, draw, drawid, group) 

mfx %>%
  ggplot(aes(x = draw, y = group, fill = after_stat(x > 1/4))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), fatten_point = 1.8,
               point_interval = "median_hdi", normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("gray80", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 3), 
                     expand = c(0, 0.02)) +
  facet_grid(. ~ gender, scales = "free_x") + 
  labs(title = "marginal effect of gender on HRF group",
       x = "probability",
       y = "HRF group") +
  theme(
    legend.position = "none",
    plot.margin = margin(5.5, 50, 5.5, 5.5, "pt"),
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(8, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = element_text(size = 32),
    text = element_text(size = 36, face = "plain", color = "black"),
    axis.text = element_text(color = "black", size = 36),
    axis.text.y = element_text(margin = margin(r = 20)),
    plot.title.position = "plot",
    plot.title = element_text(size = 46),
    plot.subtitle = element_text(size = 40),
    plot.caption = element_text(size = 40)) -> p

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 18, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Summary of the distribution

# A tibble: 8 × 8
# Groups:   gender, group [8]
  gender group Parameter Median    CI CI_low CI_high    pd
  <fct>  <fct> <chr>      <dbl> <dbl>  <dbl>   <dbl> <dbl>
1 F      1     Posterior 0.160   0.95 0.154    0.168     1
2 F      2     Posterior 0.348   0.95 0.339    0.357     1
3 F      3     Posterior 0.289   0.95 0.280    0.297     1
4 F      4     Posterior 0.203   0.95 0.195    0.210     1
5 M      1     Posterior 0.219   0.95 0.212    0.227     1
6 M      2     Posterior 0.419   0.95 0.408    0.428     1
7 M      3     Posterior 0.0959  0.95 0.0898   0.102     1
8 M      4     Posterior 0.266   0.95 0.256    0.274     1

mfx <- avg_predictions(fit1, by = "type_dur") %>% 
  posterior_draws() %>% 
  select(type_dur, draw, drawid, group) 

mfx %>%
  ggplot(aes(x = draw, y = group, fill = after_stat(x > 1/4))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), fatten_point = 1.8,
               point_interval = "median_hdi", normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("gray80", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 3), 
                     expand = c(0, 0.02)) +
  facet_grid(. ~ type_dur, scales = "free_x") + 
  labs(title = "marginal effect of IED type on HRF group",
       x = "probability",
       y = "HRF group") +
  theme(
    legend.position = "none",
    plot.margin = margin(5.5, 50, 5.5, 5.5, "pt"),
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(8, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = element_text(size = 32),
    text = element_text(size = 36, face = "plain", color = "black"),
    axis.text = element_text(color = "black", size = 36),
    axis.text.y = element_text(margin = margin(r = 20)),
    plot.title.position = "plot",
    plot.title = element_text(size = 46),
    plot.subtitle = element_text(size = 40),
    plot.caption = element_text(size = 40)) -> p

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 18, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Summary of the distribution

# A tibble: 8 × 8
# Groups:   type_dur, group [8]
  type_dur     group Parameter Median    CI CI_low CI_high    pd
  <chr>        <fct> <chr>      <dbl> <dbl>  <dbl>   <dbl> <dbl>
1 long-lasting 1     Posterior  0.184  0.95  0.178   0.189     1
2 long-lasting 2     Posterior  0.378  0.95  0.371   0.385     1
3 long-lasting 3     Posterior  0.199  0.95  0.193   0.204     1
4 long-lasting 4     Posterior  0.240  0.95  0.233   0.246     1
5 spike        1     Posterior  0.216  0.95  0.200   0.232     1
6 spike        2     Posterior  0.401  0.95  0.381   0.419     1
7 spike        3     Posterior  0.200  0.95  0.186   0.215     1
8 spike        4     Posterior  0.182  0.95  0.168   0.197     1

Within pathological cavity

HRF features

Peak amplitude

mean_posPeakAmp = mean(df$posPeakAmp)
sd_posPeakAmp = sd(df$posPeakAmp)

fit2 <- brm(
  data = df %>% mutate(posPeakAmp = (posPeakAmp - mean_posPeakAmp)/sd_posPeakAmp),
  family = gaussian(),
  formula = posPeakAmp ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakAmp_pathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakAmp_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_posPeakAmp) + mean_posPeakAmp) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "a | marginal effect of pathology on HRF peak amplitude",
       x = "BOLD % change",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Time to peak

mean_posPeakTime = mean(df$posPeakTime)
sd_posPeakTime = sd(df$posPeakTime)

fit2 <- brm(
  data = df %>% mutate(posPeakTime = (posPeakTime - mean_posPeakTime)/sd_posPeakTime),
  family = gaussian(),
  formula = posPeakTime ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakTime_pathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakTime_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_posPeakTime) + mean_posPeakTime) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "c | marginal effect of pathology on HRF peak time",
       x = "Time(s)",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Peak width

mean_posPeakFWHM = mean(df$posPeakFWHM)
sd_posPeakFWHM = sd(df$posPeakFWHM)

fit2 <- brm(
  data = df %>% mutate(posPeakFWHM = (posPeakFWHM - mean_posPeakFWHM)/sd_posPeakFWHM),
  family = gaussian(),
  formula = posPeakFWHM ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_posPeakFWHM_pathology.rds")

fit2 <- readRDS("fit_brms_model_posPeakFWHM_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_posPeakFWHM) + mean_posPeakFWHM) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))


mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "e | marginal effect of pathology on HRF peak FWHM",
       x = "FWHM(s)",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Undershoot amplitude

mean_negPeakAmp = mean(df$negPeakAmp)
sd_negPeakAmp = sd(df$negPeakAmp)

fit2 <- brm(
  data = df %>% mutate(negPeakAmp = (negPeakAmp - mean_negPeakAmp)/sd_negPeakAmp),
  family = gaussian(),
  formula = negPeakAmp ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakAmp_pathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakAmp_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_negPeakAmp) + mean_negPeakAmp) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "b | marginal effect of pathology on HRF undershoot amplitude",
       x = "BOLD % change",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Time to undershoot

mean_negPeakTime = mean(df$negPeakTime)
sd_negPeakTime = sd(df$negPeakTime)

fit2 <- brm(
  data = df %>% mutate(negPeakTime = (negPeakTime - mean_negPeakTime)/sd_negPeakTime),
  family = gaussian(),
  formula = negPeakTime ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakTime_pathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakTime_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_negPeakTime) + mean_negPeakTime) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "d | marginal effect of pathology on HRF undershoot time",
       x = "Time(s)",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Undershoot width

mean_negPeakFWHM = mean(df$negPeakFWHM)
sd_negPeakFWHM = sd(df$negPeakFWHM)

fit2 <- brm(
  data = df %>% mutate(negPeakFWHM = (negPeakFWHM - mean_negPeakFWHM)/sd_negPeakFWHM),
  family = gaussian(),
  formula = negPeakFWHM ~ age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
  prior = c(
    prior(normal(0, 1), class = b),
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 0.5), class = sd),
    prior(normal(0, 1), class = sigma)
  ),
  iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221
)

saveRDS(fit2, "fit_brms_model_negPeakFWHM_pathology.rds")

fit2 <- readRDS("fit_brms_model_negPeakFWHM_pathology.rds")

mfx <- avg_predictions(fit2, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid) %>% 
  mutate(draw = (draw * sd_negPeakFWHM) + mean_negPeakFWHM) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

mfx %>%
  ggplot(aes(x = draw, y = pathology, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("skyblue", "skyblue")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
  labs(subtitle = "f | marginal effect of pathology on HRF undershoot FWHM",
       x = "FWHM(s)",
       y = " ") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 26),
    axis.text.x = element_text(size = 24),   
    axis.text.y = ggtext::element_markdown(size = 24, margin = margin(r = 0))
  ) -> p_pathology

p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 8, units = "in", res = figdpi)
suppressWarnings(print(p_pathology))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

HRF shape

fit_brms4 <- brm(
 data = df, family = categorical(link = "logit", refcat = "1"),
 formula = hrf_group ~ 0 + age + gender + type_dur + sign + (1 | pathology) + (1 | lobe) + (1 | fs_label),
 prior = c(
   prior(normal(0, 3), class = b, dpar = mu2),
   prior(normal(0, 3), class = b, dpar = mu3),
   prior(normal(0, 3), class = b, dpar = mu4),
   prior(exponential(1), class = sd, dpar = mu2),
   prior(exponential(1), class = sd, dpar = mu3),
   prior(exponential(1), class = sd, dpar = mu4)
 ),
 control = list(adapt_delta = 0.9),
 iter = 2500, warmup = 1000, chains = 4, cores = 4, seed = 19860221)


saveRDS(fit_brms4, "fit_brms4_4groups.rds")

fit_brms4 <- readRDS("fit_brms4_4groups.rds")

fit <- fit_brms4

mfx <- avg_predictions(fit, by = "pathology") %>% 
  posterior_draws() %>% 
  select(pathology, draw, drawid, group) %>% 
  mutate(pathology = factor(pathology, levels = c("HS|dual", "cortical malformation", "complex", "normal", "gliosis", "encephalitis")))

group_colors <- c(
  "1" = "#00468BFF",
  "2" = "#ED0000FF",
  "3" = "#42B540FF",
  "4" = "#0099B4FF"
)

mfx <- mfx %>%
  mutate(group = factor(group),
         group_colored = paste0(
           "<span style='color:", group_colors[as.character(group)], "'>",
           as.character(group), "</span>"
         ))

mfx$group_colored <- factor(mfx$group_colored, levels = unique(mfx$group_colored))

mfx %>%
  ggplot(aes(x = draw, y = group_colored, fill = after_stat(x > 1/4))) +
  stat_halfeye(.width = c(.95, .5), interval_size_range = 1.5 * c(0.6, 1.4), 
               fatten_point = 1.8, point_interval = "median_hdi", 
               normalize = "groups", height = 0.8) +
  scale_fill_manual(values = c("gray80", "skyblue")) +
  geom_vline(xintercept = 0.25, linetype = "dashed") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 4)) +
  facet_wrap(~pathology, nrow = 2, ncol = 3, scales = "free_x") +  
  labs(subtitle = "c | marginal effect of pathology on HRF group",
       x = "Probability",
       y = "HRF group") +
  theme(
    plot.title.position = "plot",
    plot.subtitle = element_text(size = 32),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    panel.spacing = unit(2, "lines"),
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = ggtext::element_markdown(size = 20, hjust = 0),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 24),
    axis.text.x = element_text(size = 20),   
    axis.text.y = ggtext::element_markdown(size = 20, margin = margin(r = 20))
  ) -> p_pathology

Summary of the distribution

# A tibble: 24 × 8
# Groups:   pathology, group [24]
   pathology             group Parameter Median    CI  CI_low CI_high    pd
   <fct>                 <fct> <chr>      <dbl> <dbl>   <dbl>   <dbl> <dbl>
 1 HS|dual               1     Posterior 0.362   0.95 0.328    0.395      1
 2 HS|dual               2     Posterior 0.270   0.95 0.232    0.308      1
 3 HS|dual               3     Posterior 0.0803  0.95 0.0580   0.103      1
 4 HS|dual               4     Posterior 0.287   0.95 0.254    0.320      1
 5 cortical malformation 1     Posterior 0.0682  0.95 0.0537   0.0840     1
 6 cortical malformation 2     Posterior 0.622   0.95 0.595    0.651      1
 7 cortical malformation 3     Posterior 0.126   0.95 0.108    0.146      1
 8 cortical malformation 4     Posterior 0.183   0.95 0.159    0.206      1
 9 complex               1     Posterior 0.103   0.95 0.00831  0.276      1
10 complex               2     Posterior 0.529   0.95 0.264    0.785      1
# ℹ 14 more rows

mfx <- avg_predictions(fit, by = "fs_label")


roi_temp <- mfx %>% group_by(fs_label) %>%
  slice_max(order_by = estimate, n = 1, with_ties = FALSE) %>%
  ungroup() %>% 
  select(group, fs_label, estimate) %>% 
  mutate(fs_label = as.numeric(as.character(fs_label))) %>% 
  left_join(fs_label %>% select(fs_label, FS_labelname), by = "fs_label") %>% 
  mutate(FS_labelname = str_remove(FS_labelname, "ctx_")) %>% 
  rename(region = FS_labelname)


colorbar_col <- adjustcolor(
  c("#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f"),
  alpha.f = 1)

atlas_name <- desterieux
atlas_name$data <- atlas_name$data %>% rename(roi_name = region) %>% rename(region = label)

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  labs(subtitle = "b | posterior probability of selected HRF per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 24),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> proi1

roi_temp %>% 
  ggseg(atlas = atlas_name, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  labs(subtitle = "a | most probable HRF group per ROI") + 
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 26),             
        legend.title = element_text(size = 26),   
        axis.title = element_blank(),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 32),
        axis.text = element_text(size = 26, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF groups",
    title.position = "top",
    title.hjust = 0.5,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> proi2


subcor <- ggseg::aseg
subcor$data <- subcor$data %>% 
  filter(str_detect(label,
                    'Thalamus|Caudate|Putamen|Pallidum|Hippocampus|Amygdala|VentralDC')) %>% 
  mutate(side = "left             right")

subcormap <- subcor
subcormap$data <- subcormap$data %>% rename(roi = region) %>% rename(region = label)


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = estimate), 
        color = "black", position = "dispersed") +
  scale_fill_gradientn(limits = c(0.25, 0.9), 
                       colours = colorbar_col, 
                       na.value = "grey50", 
                       breaks = pretty_breaks(5)) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_colourbar(title = "higest probility",
                                barwidth = 25, barheight = 1, 
                                direction = "horizontal", 
                                title.position = "left", 
                                title.vjust = 1,
                                label.position = "bottom")) -> psc1


roi_temp %>% 
  ggseg(atlas = subcormap, mapping = aes(fill = group), 
        color = "black", position = "dispersed") +
  scale_fill_manual(values = c("1" = "#00468BFF", "2" = "#ED0000FF", "3" = "#42B540FF", "4" = "#0099B4FF",                        "5" = "#925E9FFF", "6" = "#FDAF91FF", "7" = "#AD002AFF", "8" = "#ADB6B6FF", "9" = "#1B1919FF")) +
  theme(legend.position = "bottom",
        legend.box = "horizontal",
        legend.text = element_text(size = 20), 
        axis.title = element_blank(),
        axis.text = element_text(size = 20, family = "", color = "black"),
        text = element_text(size = 20, face = "plain", family = "")) +
  guides(fill = guide_legend(
    title = "HRF group",
    title.position = "left",
    title.vjust = 1,
    label.position = "bottom",
    direction = "horizontal",
    nrow = 1,
    byrow = TRUE,
    keywidth = unit(2.5, "cm"),
    keyheight = unit(0.4, "cm"))) -> psc2


proisc2a <- ggarrange(proi2, psc2, nrow = 1, ncol = 2, common.legend = T, 
                    legend="bottom", widths = c(6, 1))


proisc2 <- ggarrange(proisc2a, proisc2b, p_pathology, nrow = 3, ncol = 1, 
                     heights = c(0.9, 1, 1.5), align = "hv")


p_temp_file <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file, width = 16, height = 2*4.5+8, units = "in", res = figdpi)
suppressWarnings(print(proisc2))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

Model comparsion

Age, gender, IED type and HRF sign complete pooling

fit0 <- mblogit(
  formula = hrf_group ~ age + gender + type_dur + sign,
  data    = df,
  catCov  = "free",      
  method  = "PQL",         
  estimator = "ML"
)

Model fitting

Iteration 1 - deviance = 3504.293 - criterion = 0.3852283
Iteration 2 - deviance = 3479.118 - criterion = 0.007235657
Iteration 3 - deviance = 3478.492 - criterion = 0.0001799623
Iteration 4 - deviance = 3478.491 - criterion = 2.988217e-07
Iteration 5 - deviance = 3478.491 - criterion = 1.498659e-12
converged

Age, gender, IED type and HRF sign complete pooling. Parcel and lobe two level hierarchic partial pooling.

fit1 <- mblogit(
  formula = hrf_group ~ age + gender + type_dur + sign,
  data    = df,
  random  = list(~1 | lobe,
                 ~1 | fs_label),
  catCov  = "free",      
  method  = "PQL",         
  estimator = "ML"
)

Model fitting

Iteration 1 - deviance = 3161.453 - criterion = 0.8412412
Iteration 2 - deviance = 3103.597 - criterion = 0.0532227
Iteration 3 - deviance = 3098.86 - criterion = 0.005489007
Iteration 4 - deviance = 3096.22 - criterion = 7.67433e-05
Iteration 5 - deviance = 3096.07 - criterion = 1.952338e-08
Iteration 6 - deviance = 3096.069 - criterion = 1.789622e-15
converged

Age, gender, IED type and HRF sign complete pooling. Pathology hierarchic partial pooling.

fit2 <- mblogit(
  formula = hrf_group ~ age + gender + type_dur + sign,
  data    = df,
  random  = list(~1 | pathology),
  catCov  = "free",      
  method  = "PQL",         
  estimator = "ML"
)

Model fitting

Iteration 1 - deviance = 3315.923 - criterion = 0.8989671
Iteration 2 - deviance = 3284.55 - criterion = 0.05741227
Iteration 3 - deviance = 3282.867 - criterion = 0.003948137
Iteration 4 - deviance = 3282.538 - criterion = 0.0002067492
Iteration 5 - deviance = 3282.492 - criterion = 1.004058e-07
Iteration 6 - deviance = 3282.49 - criterion = 2.315406e-14
converged

Age, gender and IED type complete pooling. Parcel and lobe two level hierarchic partial pooling. Pathology hierarchic partial pooling.

fit3 <- mblogit(
  formula = hrf_group ~ age + gender + type_dur,
  data    = df,
  random  = list(~1 | pathology,
                 ~1 | lobe,
                 ~1 | fs_label),
  catCov  = "free",      
  method  = "PQL",         
  estimator = "ML"
)

Model fitting

Iteration 1 - deviance = 3047.633 - criterion = 0.8392789
Iteration 2 - deviance = 2973.152 - criterion = 0.06376108
Iteration 3 - deviance = 2964.011 - criterion = 0.01196605
Iteration 4 - deviance = 2960.176 - criterion = 0.0006057058
Iteration 5 - deviance = 2959.683 - criterion = 1.897951e-06
Iteration 6 - deviance = 2959.665 - criterion = 2.148408e-11
converged

Age, gender, IED type and HRF sign complete pooling. Parcel and lobe two level hierarchic partial pooling. Pathology hierarchic partial pooling.

fit4 <- mblogit(
  formula = hrf_group ~ age + gender + type_dur + sign,
  data    = df,
  random  = list(~1 | pathology,
                 ~1 | lobe,
                 ~1 | fs_label),
  catCov  = "free",      
  method  = "PQL",         
  estimator = "ML"
)

Model fitting

Iteration 1 - deviance = 3039.761 - criterion = 0.8405424
Iteration 2 - deviance = 2965.393 - criterion = 0.06316721
Iteration 3 - deviance = 2956.404 - criterion = 0.00895744
Iteration 4 - deviance = 2953.072 - criterion = 0.0003377023
Iteration 5 - deviance = 2952.711 - criterion = 7.758474e-07
Iteration 6 - deviance = 2952.7 - criterion = 5.216892e-12
converged

Summary of the performance

# Comparison of Model Performance Indices

Name  |                                                                  Model
------------------------------------------------------------------------------
rank1 |        hrf ~ age + gender + type + (1 | lobe/parcel) + (1 | pathology)
rank2 | hrf ~ age + gender + type + sign + (1 | lobe/parcel) + (1 | pathology)
rank3 |                   hrf ~ age + gender + type + sign + (1 | lobe/parcel)
rank4 |                     hrf ~ age + gender + type + sign + (1 | pathology)
rank5 |                                       hrf ~ age + gender + type + sign

Name  | Nagelkerke.s.R2 |  RMSE | Sigma | AIC weights | BIC weights | Performance-Score
---------------------------------------------------------------------------------------
rank1 |           0.619 | 0.340 | 1.385 |       0.382 |       0.999 |            92.08%
rank2 |           0.621 | 0.340 | 1.385 |       0.618 |    5.31e-04 |            79.92%
rank3 |           0.578 | 0.348 | 1.418 |    1.84e-29 |    1.47e-25 |            46.17%
rank4 |           0.516 | 0.371 | 1.460 |    2.45e-67 |    1.82e-56 |            21.94%
rank5 |           0.441 | 0.386 | 1.503 |   2.72e-107 |    1.88e-89 |             0.00%

Classifications

Classifying pathology

df_est <- df %>%
  filter(pathology != "complex") %>%
  mutate(pathology = factor(pathology))

set.seed(19860221)

# define nested resamples (5-fold outer, 5-fold inner, stratified)

nested_splits <- nested_cv(
  df_est,
  outside = vfold_cv(v = 5, strata = pathology),
  inside  = vfold_cv(v = 5, strata = pathology)
)


# recipe

rec <- recipe(pathology ~ posPeakAmp + posPeakTime + posPeakFWHM +
                          negPeakAmp + negPeakTime + negPeakFWHM,
              data = df_est) %>%
  step_upsample(pathology, over_ratio = 0.5) %>%
  step_dummy(all_nominal_predictors())

Classifiers were evaluated on the pre-specified 20% stratified hold-out set using the same inverse-frequency case weights derived from the training set. We extracted predicted and true labels and showed a confusion-matrix first.

We quantified agreement and discrimination with five complementary measures:

• Cohen’s κ: chance-corrected agreement across all five classes,

• Matthews Correlation Coefficient (MCC): a balanced summary of true/false positives and negatives,

• Balanced Accuracy: the average of per-class recall, which prevents dominant classes from inflating our score,

• Macro-F₁: the unweighted mean of each class’s F₁ score, and

• Weighted-F₁: the prevalence-weighted mean of F₁ to reflect the observed class distribution.

To access class-specific performance, we computed one-versus-all versions of Cohen’s κ, MCC, balanced accuracy, and F₁ for each pathology. This breakdown highlights which pathologies the model separates reliably.

# workflow

rf_spec <- rand_forest(
  trees = 500,
  mtry  = tune(),
  min_n = tune()
) %>%
  set_engine(
    "ranger",
    probability = T,
    importance = "impurity"
  ) %>%
  set_mode("classification")


rf_wf <- workflow() %>%
  add_recipe(rec) %>%
  add_model(rf_spec)


# parameter grid

rf_params <- parameters(rf_spec) %>%
  update(
    mtry  = mtry(c(1, sum(rec$var_info$role == "predictor"))),
    min_n = min_n(c(1, round(nrow(df_est)*0.8*0.8*0.01)))
  )

rf_grid <- grid_random(rf_params, size = 100)


# inner-loop tuning, outer-loop fitting

tuned_nested <- nested_splits %>%
  mutate(
    tune_res = map(inner_resamples, ~ tune_grid(
      rf_wf,
      resamples = .x,
      grid = rf_grid,
      metrics = metric_set(f_meas),
      control = control_resamples(save_pred = T)
    )),
    best_params = map(tune_res,  ~ select_best(.x, metric = "f_meas")),
    final_wf = map(best_params, ~ finalize_workflow(rf_wf, .x)),
    outer_fit = map2(final_wf, splits, ~ last_fit(.x, .y)),
    preds = map(outer_fit, ~ collect_predictions(.x))
  )


# pool all outer‐holdout predictions

all_preds <- tuned_nested %>%
  select(preds) %>%
  unnest(preds)

# confusion matrix

cm <- conf_mat(all_preds, truth = pathology, estimate = .pred_class)

mat_counts <- as.matrix(cm$table)
col_sums <- colSums(mat_counts)
mat_norm  <- sweep(mat_counts, 2, col_sums, FUN = "/")

df_counts <- as.data.frame(as.table(mat_counts))
names(df_counts) <- c("Prediction", "Truth", "Count")

df_norm <- as.data.frame(as.table(mat_norm))
names(df_norm) <- c("Prediction", "Truth", "Freq")

df_plot <- left_join(df_norm, df_counts, by = c("Prediction", "Truth")) %>%
   mutate(Prediction = as.character(Prediction),
          Truth = as.character(Truth))

desired_order <- c(
  "HS|dual",
  "cortical malformation",
  "normal",
  "gliosis",
  "encephalitis"
)
df_plot$Prediction <- factor(df_plot$Prediction, levels = desired_order)
df_plot$Truth <- factor(df_plot$Truth, levels = desired_order)


ggplot(df_plot, aes(x = Truth, y = Prediction, fill = Freq)) +
  geom_tile(color = "white") +
  scale_fill_material("red") +
  scale_y_discrete(limits = rev(desired_order)) +
  geom_text(aes(label = Count), color = "black", size = 6) +
  labs(x = "Actual", y = "Classified") +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    legend.position = "none",
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    axis.ticks = element_blank(),
    axis.title.y = element_text(size = 20),
    axis.title.x = element_text(size = 20),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size = 20)
  ) +
  coord_fixed() -> p


p_temp_file1 <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file1, width = 5.33, height = 5.33, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())


# metrics

# overall

fold_metrics <- map2_dfr(
  tuned_nested$preds,
  seq_along(tuned_nested$preds),
  function(preds, fold) {
    tibble(
      fold = fold,
      kappa = kap_vec(preds$pathology, preds$.pred_class),
      mcc = mcc_vec(preds$pathology, preds$.pred_class),
      precision = precision_vec(preds$pathology, preds$.pred_class, estimator = "macro"),
      recall  = recall_vec(preds$pathology, preds$.pred_class, estimator = "macro"),
      F1  = f_meas_vec(preds$pathology, preds$.pred_class, estimator = "macro")
    )
  }
)

fold_metrics %>%
  summarise(across(-fold, list(mean = mean, sd = sd), .names = "{.col}_{.fn}")) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("metric", "stat"),
    names_pattern = "(.*)_(mean|sd)"
  ) %>%
  pivot_wider(
    names_from = stat,
    values_from  = value
  ) %>%
  mutate(mean_sd = sprintf("%.2f ± %.2f", mean, sd)) %>%
  select(metric, mean_sd) %>%
  select(metric, mean_sd) %>%
  pivot_wider(
    names_from = metric,
    values_from = mean_sd
  ) %>%
  kbl(format = "html", caption = "overall metrics (mean ± sd over outer folds)") %>%
  kable_paper(full_width = TRUE)

overall metrics (mean ± sd over outer folds)

kappa	mcc	precision	recall	F1
0.74 ± 0.03	0.74 ± 0.03	0.79 ± 0.04	0.79 ± 0.05	0.79 ± 0.04

fold_importance <- tuned_nested %>%
  mutate(
    vi = map(outer_fit, ~ {
      wf_fitted  <- extract_workflow(.x)
      ranger_obj <- pull_workflow_fit(wf_fitted)$fit
      tibble(
        variable = names(ranger_obj$variable.importance),
        importance = ranger_obj$variable.importance
      )
    })
  ) %>%
  select(vi) %>%
  unnest(vi, .id = "fold")


mean_importance <- fold_importance %>%
  group_by(variable) %>%
  summarise(
    mean_imp = mean(importance),
    sd_imp   = sd(importance),
    .groups  = "drop"
  ) %>%
  arrange(desc(mean_imp))

mean_importance %>%
  mutate(variable = recode(variable,
    posPeakAmp = "Amplitude (peak)",
    posPeakTime = "Time (peak)",
    posPeakFWHM = "FWHM (peak)",
    negPeakAmp = "Amplitude (undershoot)",
    negPeakTime = "Time (undershoot)",
    negPeakFWHM = "FWHM (undershoot)"
  )) %>%
  slice_max(order_by = mean_imp, n = 20) %>%
  mutate(variable = factor(variable, levels = rev(variable))) %>%
  ggplot(aes(x = variable, y = mean_imp)) +
  geom_col(fill = "skyblue") +
  geom_text(
      aes(x = variable, y = mean_imp/2, label = variable),
      color = "black",
      size = 6
    ) +
  geom_errorbar(aes(ymin = mean_imp - sd_imp,
                      ymax = mean_imp + sd_imp),
                  width = 0.2) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
  coord_flip() +
  labs(
      title = "Variable importance (mean ± sd)",
      x = NULL,
      y = "Gini importance"
    ) +
  theme_minimal() +
  theme(
      axis.text.y = element_blank(),
      panel.background = element_blank(),
      panel.grid = element_blank(),
      legend.position = "none",
      plot.title = element_text(size = 22, color = "black"),
      axis.title.x = element_text(size = 22, color = "black"),
      axis.text.x = element_text(size = 20, color = "black")
    ) -> p


p_temp_file2 <- tempfile(tmpdir = "./tmp", fileext = '.png')
agg_png(p_temp_file2, width = 6, height = 6, units = "in", res = figdpi)
suppressWarnings(print(p))
invisible(dev.off())

Classifying lobe

df_est <- df %>%
  mutate(
    lobe = str_remove(lobe, "^(left |right )")
  ) %>%
  mutate(lobe = factor(lobe))

set.seed(19860221)

# define nested resamples (5-fold outer, 5-fold inner, stratified)

nested_splits <- nested_cv(
  df_est,
  outside = vfold_cv(v = 5, strata = lobe),
  inside  = vfold_cv(v = 5, strata = lobe)
)


# recipe

rec <- recipe(lobe ~ posPeakAmp + posPeakTime + posPeakFWHM +
                negPeakAmp + negPeakTime + negPeakFWHM,
              data = df_est) %>%
  step_upsample(lobe, over_ratio = 0.5) %>%
  step_dummy(all_nominal_predictors())

# workflow

rf_spec <- rand_forest(
  trees = 500,
  mtry = tune(),
  min_n = tune()
) %>%
  set_engine(
    "ranger",
    probability = T,
    importance = "impurity"
  ) %>%
  set_mode("classification")


rf_wf <- workflow() %>%
  add_recipe(rec) %>%
  add_model(rf_spec)


# parameter grid

rf_params <- parameters(rf_spec) %>%
  update(
    mtry = mtry(c(1, sum(rec$var_info$role == "predictor"))),
    min_n = min_n(c(1, round(nrow(df_est)*0.8*0.8*0.01)))
  )

rf_grid <- grid_random(rf_params, size = 100)


# inner-loop tuning, outer-loop fitting

tuned_nested <- nested_splits %>%
  mutate(
    tune_res = map(inner_resamples, ~ tune_grid(
      rf_wf,
      resamples = .x,
      grid = rf_grid,
      metrics = metric_set(f_meas),
      control = control_resamples(save_pred = T)
    )),
    best_params = map(tune_res,  ~ select_best(.x, metric = "f_meas")),
    final_wf = map(best_params, ~ finalize_workflow(rf_wf, .x)),
    outer_fit = map2(final_wf, splits, ~ last_fit(.x, .y)),
    preds = map(outer_fit, ~ collect_predictions(.x))
  )


# pool all outer‐holdout predictions

all_preds <- tuned_nested %>%
  select(preds) %>%
  unnest(preds)

# metrics

fold_metrics <- map2_dfr(
  tuned_nested$preds,
  seq_along(tuned_nested$preds),
  function(preds, fold) {
    tibble(
      fold = fold,
      kappa = kap_vec(preds$lobe, preds$.pred_class),
      mcc = mcc_vec(preds$lobe, preds$.pred_class),
      precision = precision_vec(preds$lobe, preds$.pred_class, estimator = "macro"),
      recall = recall_vec(preds$lobe, preds$.pred_class, estimator = "macro"),
      F1 = f_meas_vec(preds$lobe, preds$.pred_class, estimator = "macro")
    )
  }
)

fold_metrics %>%
  summarise(across(-fold, list(mean = mean, sd = sd), .names = "{.col}_{.fn}")) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("metric", "stat"),
    names_pattern = "(.*)_(mean|sd)"
  ) %>%
  pivot_wider(
    names_from = stat,
    values_from = value
  ) %>%
  mutate(mean_sd = sprintf("%.2f ± %.2f", mean, sd)) %>%
  select(metric, mean_sd) %>%
  select(metric, mean_sd) %>%
  pivot_wider(
    names_from = metric,
    values_from = mean_sd
  ) %>%
  kbl(format = "html", caption = "overall metrics (mean ± sd over outer folds)") %>%
  kable_paper(full_width = TRUE)

overall metrics (mean ± sd over outer folds)

kappa	mcc	precision	recall	F1
0.56 ± 0.05	0.56 ± 0.05	0.62 ± 0.06	0.55 ± 0.06	0.57 ± 0.05

Computing environment

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS 15.7.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Toronto
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] vip_0.4.1                 ranger_0.17.0            
 [3] xgboost_1.7.11.1          themis_1.0.3             
 [5] yardstick_1.3.2           workflowsets_1.1.1       
 [7] workflows_1.2.0           tune_1.3.0               
 [9] rsample_1.3.0             recipes_1.3.1            
[11] parsnip_1.3.2             modeldata_1.4.0          
[13] infer_1.0.8               dials_1.4.0              
[15] broom_1.0.8               tidymodels_1.3.0         
[17] ggridges_0.5.6            ggsegSchaefer_2.0.4      
[19] ggsegDesterieux_1.0.1.002 ggseg_1.6.5              
[21] mclogit_0.9.6             lme4_1.1-37              
[23] Matrix_1.7-0              neuRosim_0.2-14          
[25] deSolve_1.40              nnet_7.3-19              
[27] brainspriteR_1.0.0.9000   ggrepel_0.9.6            
[29] marginaleffects_0.27.0    ggsci_3.2.0              
[31] ggforce_0.5.0             moonBook_0.3.1           
[33] kableExtra_1.4.0          ggdag_0.2.13             
[35] dagitty_0.3-4             DT_0.33                  
[37] webr_0.1.5                tidybayes_3.0.7          
[39] rstan_2.32.7              StanHeaders_2.32.10      
[41] brms_2.22.0               Rcpp_1.0.14              
[43] ggstatsplot_0.13.1        gt_1.0.0                 
[45] readxl_1.4.5              ragg_1.4.0               
[47] scales_1.4.0              ggpubr_0.6.0             
[49] lubridate_1.9.4           forcats_1.0.0            
[51] stringr_1.5.1             dplyr_1.1.4              
[53] purrr_1.0.4               readr_2.1.5              
[55] tidyr_1.3.1               tibble_3.3.0             
[57] ggplot2_4.0.0             tidyverse_2.0.0          

loaded via a namespace (and not attached):
  [1] vroom_1.6.5             sjmisc_2.8.10           vctrs_0.6.5            
  [4] effectsize_1.0.1        digest_0.6.37           proxy_0.4-27           
  [7] bayestestR_0.16.1       correlation_0.8.7       parallelly_1.45.0      
 [10] MASS_7.3-60.2           fontLiberation_0.1.0    httpuv_1.6.16          
 [13] foreach_1.5.2           withr_3.0.2             psych_2.5.6            
 [16] xfun_0.52               survival_3.6-4          commonmark_1.9.5       
 [19] emmeans_1.11.1          parameters_0.26.0       systemfonts_1.2.3      
 [22] zoo_1.8-14              V8_6.0.4                ggdist_3.3.3           
 [25] Formula_1.2-5           datawizard_1.1.0        rematch2_2.1.2         
 [28] promises_1.3.3          httr_1.4.7              rstatix_0.7.3          
 [31] globals_0.18.0          ps_1.9.1                rstudioapi_0.17.1      
 [34] units_0.8-7             miniUI_0.1.2            generics_0.1.4         
 [37] processx_3.8.6          curl_6.4.0              polyclip_1.10-7        
 [40] xtable_1.8-4            svUnit_1.0.6            evaluate_1.0.4         
 [43] hms_1.1.3               ggseg3d_1.6.3           colorspace_2.1-1       
 [46] shinyWidgets_0.9.0      magrittr_2.0.3          lmtest_0.9-40          
 [49] later_1.4.4             lattice_0.22-6          posterior_1.6.1        
 [52] future.apply_1.20.0     lhs_1.2.0               cowplot_1.1.3          
 [55] matrixStats_1.5.0       ROSE_0.0-4              class_7.3-22           
 [58] pillar_1.10.2           nlme_3.1-164            performance_0.14.0     
 [61] iterators_1.0.14        GPfit_1.0-9             compiler_4.4.1         
 [64] stringi_1.8.7           gower_1.0.2             sf_1.0-21              
 [67] minqa_1.2.8             crayon_1.5.3            abind_1.4-8            
 [70] ggtext_0.1.2            bit_4.6.0               codetools_0.2-20       
 [73] textshaping_1.0.1       openssl_2.3.3           flextable_0.9.9        
 [76] bslib_0.9.0             QuickJSR_1.8.0          e1071_1.7-16           
 [79] paletteer_1.6.0         plotly_4.11.0           mime_0.13              
 [82] sparsevctrs_0.3.4       splines_4.4.1           markdown_2.0           
 [85] DiceDesign_1.10         cellranger_1.1.0        gridtext_0.1.5         
 [88] utf8_1.2.6              knitr_1.50              listenv_0.9.1          
 [91] checkmate_2.3.2         Rdpack_2.6.4            pkgbuild_1.4.8         
 [94] openxlsx_4.2.8          ggsignif_0.6.4          estimability_1.5.1     
 [97] callr_3.7.6             tzdb_0.5.0              svglite_2.2.1          
[100] tweenr_2.0.3            bayesplot_1.13.0        pkgconfig_2.0.3        
[103] tools_4.4.1             cachem_1.1.0            rbibutils_2.3          
[106] viridisLite_0.4.2       DBI_1.2.3               fastmap_1.2.0          
[109] rmarkdown_2.29          grid_4.4.1              sass_0.4.10            
[112] officer_0.6.10          patchwork_1.3.1         coda_0.19-4.1          
[115] insight_1.3.1           carData_3.0-5           rpart_4.1.23           
[118] farver_2.1.2            reformulas_0.4.1        tidygraph_1.3.1        
[121] yaml_2.3.10             cli_3.6.5               stats4_4.4.1           
[124] lifecycle_1.0.4         askpass_1.2.1           mvtnorm_1.3-3          
[127] lava_1.8.1              backports_1.5.0         Brobdingnag_1.2-9      
[130] timechange_0.3.0        gtable_0.3.6            arrayhelpers_1.1-0     
[133] devEMF_4.5-1            parallel_4.4.1          jsonlite_2.0.0         
[136] memisc_0.99.31.8.3      bit64_4.6.0-1           see_0.11.0             
[139] litedown_0.7            zip_2.3.3               RcppParallel_5.1.10    
[142] jquerylib_0.1.4         bridgesampling_1.1-2    loo_2.8.0              
[145] zeallot_0.2.0           distributional_0.5.0    timeDate_4041.110      
[148] lazyeval_0.2.2          shiny_1.10.0            collapse_2.1.2         
[151] htmltools_0.5.8.1       rvg_0.3.5               sjlabelled_1.2.0       
[154] glue_1.8.0              editData_0.1.8          gdtools_0.4.2          
[157] classInt_0.4-11         mnormt_2.1.1            gridExtra_2.3          
[160] boot_1.3-30             igraph_2.1.4            R6_2.6.1               
[163] vcd_1.4-13              labeling_0.4.3          ztable_0.2.3           
[166] ipred_0.9-15            nloptr_2.2.1            rstantools_2.4.0       
[169] tidyselect_1.2.1        tensorA_0.36.2.1        rrtable_0.3.0          
[172] xml2_1.3.8              inline_0.3.21           fontBitstreamVera_0.1.1
[175] car_3.1-3               future_1.58.0           statsExpressions_1.7.0 
[178] KernSmooth_2.23-24      S7_0.2.0                furrr_0.3.1            
[181] rio_1.2.3               fontquiver_0.2.1        data.table_1.17.6      
[184] htmlwidgets_1.6.4       RColorBrewer_1.1-3      rlang_1.1.6            
[187] uuid_1.2-1              hardhat_1.4.1           prodlim_2025.04.28